'''
Created on Oct 29, 2009

@author: mkiyer
'''

from veggie.gene.gene import parse_refseq_genes

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import math


def histogram(data, xlabel='x', ylabel='y', outfile=None):
    data = np.array(data)
    print 'average', np.average(data)

    fig = plt.figure()
    ax = fig.add_subplot(111)
    nbins=100
    # the histogram of the data
    n, bins, patches = ax.hist(data, bins=nbins, facecolor='green', alpha=0.5, log=False, normed=0)
    plt.axvline(x=np.average(data), linewidth=4, color='r')
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    #n, bins, patches = ax.hist(x2, bins=bins, facecolor='red', alpha=0.25)
    #l, = plt.plot(bins, mlab.normpdf(bins, 0.0, 1.0), 'r--', label='fit', linewidth=3)
    #legend([l, patches[0]], ['fit', 'hist'])
    ax.grid(True)
    if outfile is None:
        plt.show()
    else:
        plt.savefig(outfile)
        
    plt.close()
    return

if __name__ == '__main__':    
    data = []
    n = 0
    x = 0
    for g in parse_refseq_genes():
        for i in g.introns:
            length = i[1] - i[0]
            if length > 60:
                n += 1
                continue
            x += 1
            data.append(i[1] - i[0])
    histogram(data)
    print x
    print n